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Bogoliubov-de Gennes theory is used to investigate the effect of the size of a superconducting 
square on the vortex states in the quantum confinement regime. When the superconducting coher- 
ence length is comparable to the Fermi wavelength, the shape resonances of the superconducting 
order parameter have strong influence on the vortex conflguration. Several unconventional vortex 
states, including asymmetric ones, giant multi-vortex combinations, and states comprising giant 
antivortex, were found as ground states and their stability was found to be very sensitive on the 
value of /cfCo, the size of the sample W , and the magnetic flux $. By increasing the temperature 
and/or enlarging the size of the sample, quantum conflnement is suppressed and the conventional 
mesoscopic vortex states as predicted by the Ginzburg-Laudau (GL) theory are recovered. However, 
contrary to the GL results we found that the states containing symmetry-induced vortex-antivortex 
pairs are stable over the whole temperature range. It turns out that the inhomogeneous order 
parameter induced by quantum conflnement favors vortex-antivortex molecules, as well as giant 
vortices with a rich structure in the vortex core - unattainable in the GL domain. 
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I. INTRODUCTION 

Vortex states in mesoscopic superconductors have been 
extensively studied in the past two decades, both theo- 
retically and experimentally. Two main interactions have 
been found to govern vortex behavior in a mesoscopic 
systemji"^^ The first one is the vortex-vortex interaction, 
which causes vortices to form Abrikosov triangular lat- 
tices in bulk type-II superconductors. The second one is 
the interaction between vortices and sample boundaries, 
which makes vortex configurations strongly dependent on 
the size and geometry of mesoscopic samples-whose di- 
mensions are of the order of the penetration depth A or 
the coherence length i^. For example, in square meso- 
scopic samples, vortex configurations try to best match 
the C4 symmetry. When there is only one vortex in the 
sample {L = 1 state where L is the winding number 
or vorticity), the vortex always sits in the center of the 
sample in order to balance the boundary effect from all 
sides. For the L = 2 state, two vortices sit on the diag- 
onal such that the vortex-vortex separation is maximal 
in order to minimize the intervortex interaction. A giant 
vortex with L = 2 can be induced when the boundary 
confinement pushes two single vortices together, as pre- 
dicted theoretically^ and observed experimentally. '-'^d^ 
For L = 3 state, because of its incomparability with the 
four-fold symmetry, the theory predicts that the ground 
state corresponds to an anti-vortex sitting at the center 
surrounded by four vorticesplSii^ In short, the symmetry 
of the sample largely determines the vortex configura- 
tions when the size of the superconductor is reduced. 

However, the properties of nanoscale superconductors, 
whose sizes are of the order of the Fermi wavelength Xp, 
are very different from those of mesoscopic superconduc- 
tors. This is because the distance between electronic lev- 
els becomes comparable to the superconducting energy 
gap due to quantum confinement^. As a consequence. 



the number of Cooper pairs is suppressed which leads to 
quantum-size effect (QSE) )^^d^ quantum-size cascades^^ 
the shell effect^ and inhomogeneous spatial distribution 
of the order parameter .'iS The latter is the most im- 
portant for the present work because it is expected to 
strongly influence the vortex states in nanoscale super- 
conductors. 

Inhomogeneous superconductivity has been studied in 
various systems in the last decades and shows more com- 
plex behavior than homogeneous ones. It is known that 
vortices tend to migrate and get pinned in areas where 
superconductivity is suppressed!^ The reason is that it 
is more favorable energetically for a vortex to suppress 
the superconducting order parameter in a region where 
it has already been suppressed, although sometimes vor- 
tices can be pinned where the gap is large 1^ Some three- 
dimensional (3D) samples can also be treated as inhomo- 
geneous systems. ^■^"^^ For example for a 3D tip geometry, 
an asymmetric L = 1 vortex state can be the ground state 
because the thick region prevents the vortex from pen- 
etrating it.— In multi-layered superconductors, vortices 
enter first and reside favorably in the weak layers. Then, 
vortices will penetrate into the strong layers only after 
weak layers become saturated and various vortex clus- 
ters and asymmetric vortex states are induced.^ Also, 
the fabrication of anti-dots in superconductors results in 
spatially varying superconducting energy gap with a bar- 
rier at the interfaces. In these systems, the combination 
of the giant vortex, multi-vortex and anti-vortex states 
can be found as ground state, which depends strongly on 
detailed geometry of the antidots*^"— 

Recently, we investigated'^'' the vortex states in nano- 
scale superconducting squares using the Bogoliubov-de 
Gennes (BdG) method. We found unconventional vortex 
states in the quantum limit due to shape-induced res- 
onances in the inhomogeneous Cooper-pair condensate. 
Vortex-antivortex structures, asymmetric vortex states 
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and vortex clusters were found as ground states over a 
wide range of parameters. They are distinct from previ- 
ous results obtained in mesoscopic superconductors using 
the GL theory. However, there are still several aspects 
that remained unclear. For example, how does the size 
of the sample affect the vortex states in the nanoscale 
superconductors? Under which conditions, does one re- 
cover the conventional GL results? Why are the antivor- 
tex states more stable in the nanoscale limit while giant 
vortex states are unfavorable? How the vortex states 
change if temperature is increased? 

In order to answer these questions, we extend our work 
in this paper on the vortex states in nanoscale supercon- 
ducting squares and provide a systematic study. Vortex 
states for different sample sizes, kp^o parameters and 
temperatures T are investigated and the stability of the 
symmetry- induced vortex/anti- vortex molecules is dis- 
cussed. More unconventional states, very different from 
the ones obtained within GL theory, are found. This 
study is therefore fully complementary to what is known 
for vortex matter in superconductors. 

The paper is organized as follows. In Sec. II, we intro- 
duce the theoretical approach, i.e. the BdG approach for 
a square geometry. In Sec. HI, we present the inhomoge- 
neous superconducting state in the absence of the mag- 
netic field in order to better understand the QSE in nano- 
scale superconductors. In Sec. IV, the ground states and 
metastable states are studied at zero temperature and 
on the sample size dependence of the vortex states is dis- 
cussed. In Sec. V, the finite temperature ground states 
are studied. In Sec. VI, we discuss the generation of 
vortex/anti- vortex molecules and study the structure of 
the vortex core. Finally, we summarize our findings in 
Sec. VH. 



II. THEORETICAL FORMALISM 

In the presence of a magnetic field, the Bogoliubov-de 
Gcnncs (BdG) equations read 

A(r)\ fun{r)\ _ (u^{r) 
,A(r)* -H*J U(^iy Un(f) 

where Un{vn) are electron- (hole-)like quasi-particle eigen- 
wave functions and En are the quasi-particle eigen- 
energies. The single-electron Hamiltonian reads = 

- —Y - Ep with Ep being the Fermi energy and 



(1) 



A the vector potential (we consider a gauge such that 
V- A = 0). Furthermore, we take the electron band-mass 
to be isotropic (i.e., nix = my = ruz = m). The pair 
potential is determined self-consistently from the eigen- 
wave functions and eigen-energies: 



A(f) 



Un{r)vl{f)[l-2fn] 



(2) 



where g is the coupling constant, Ec is the cutoff energy, 
and /n = [1 + exp(i5„/A:BT)]~^ is the Fermi distribu- 
tion function, where T is the temperature. We consider 
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I 


11 


III 


Ef/Fo 


4 


9 


25 


Ec/Fo 


30/7r 


30/tt 


50/tv 


Ao/Fo 


1.245 


1.85 


3.14 




2 


3 


5 



TABLE I. The parameters of the considered sample. Coupling 
constant g = 0.4343 and Ep is measured from the bottom of 
the quadratic band. 



a two-dimensional problem and assume a circular Fermi 
surface. The confinement imposes Dirichlet boundary 
conditions (i.e. 'u„(r) — 0, w„(r) — 0, f G dS) such 
that the order parameter vanishes at the surface. In an 
extreme type-II superconductor (and/or very thin sam- 
ple), it is reasonable to neglect the contribution of the 
supercurrent to the total magnetic field. 

In order to solve Eqs. ([TJ and ^ numerically, we 
expanded the wave functions u„ and u„ as 



u„ir) 



where the basis set 



ifiir iv(x,v) ~ — Sin sm — — 



(4) 



is the cigcn-states of the Hamiltonian in the absence 
of the magnetic field where W is the size of the square 
sample. The eigen-energies of such states are Tj^j^ = 

^iw^'^^^'x^^y)^EF- Through comparison with results 
obtained by using the finite difference method, we found 
that the results can be converged when we include the 
states with energies as large as ^/^o + 5 times the cutoff 
energy E^ i.e. Tj^j^ e (*/$o + 5) x [-Ec^Ec] where $o 
is the flux quantum and $ the flux through our sample. 
The free energy^i^ of the system is given by 

F^Y. P^"/" + ^bT[U In /n + (1 - /„) ln(l - /„)]] 
dr -2^£;„|t;„ 



2A(r) Vu>„[l-2/„ 



f^u>„(l - 2/„) 



(5) 



where the spatial dependence of u„ and w„ is implicit. 

In our numerical investigations, we restrict ourselves to 
the three materials (or samples) with parameters given 
in Table HI For convenience, we measure the distance in 
units of the bulk coherence length at zero temperature 
^0 and the energy in units of Fq = /i^/2m^Q. Here, = 
hvp /it Aq where vp is the Fermi velocity and Aq is the 
bulk value of the order parameter at zero temperature. 
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III. SPONTANEOUS INHOMOGENEOUS 
SUPERCONDUCTIVITY INDUCED BY 
QUANTUM CONFINEMENT 



First, we study the spatial distribution of the order 
parameter as a function of sample size in the absence 
of the magnetic field since the vortex configurations will 
be strongly affected by that distribution once the mag- 
netic field is applied. In Fig. [T] contour plots of the order 
parameter A(x,y) and the corresponding diagonal pro- 
files are shown for sample I-III with different sizes. As 
expected, all the order parameters are fourfold symmet- 
ric and show Friedel-like oscillations in space which re- 
sult in four well defined peaks at each corner. Hence, 
the superconductivity is inhomogeneous. For example. 
Figs. [T] (a)-(c) show the results for sample I with size 
W/^Q = 5,8,10, respectively. For W/£_o = 5, there are 
three oscillations in the order parameter along the diag- 
onal and the resonant central peak prevents vortex from 
sitting here. However, the profile of the order parameter 
can be changed dramatically when the size of the sample 
changes. For W/£,o = 8 [shown in Figs. [Ijb)], we found 
that the central peak disappears and a relatively fiatter 
area generates in the center. When W/£,o is increased to 
10, the flat area enlarges and the Friedel-like oscillations 
can be neglected at center when we compared it with 
oscillations near boundary. It indicates that the oscil- 
lations of the order parameter result from the quantum 
confinement effect (or boundary effect). 



In order to study the quantum size effect on the order 
parameter, we show the amplitude of the order param- 
eter in the center of the sample |Ac| and the spatially 
averaged value over the whole sample |A„j| for the sam- 
ple I with sizes W/£,o = 5 - 20 in Fig. H As seen, |Ac| 
changes dramatically with W/(,o increasing and converge 
to A/Fq = 1.245 when W/^q > 15. At the same time, the 
I Am I increases gradually with W/S,o increasing. In prin- 
ciple, both of the parameters will converge as — oo 
where the boundary effect can be totally neglected. Since 
I Ac I and \Am\ show strong quantum size effect between 
W^/Co = 5 and 10, we limit ourselves in following sections 
to study the samples for these particular size . 



The profile of A is also strongly affected by fci?^o- 
Figs, md) and (e) show results for sample II with size 
^/Co = 5 and 10 and (f) is for sample III with size 
^/Co = 5. Comparing to the sample I with the same 
size, we found that the wave number and the amplitude 
of the oscillations along diagonal are larger. It indicates 
that the superconducting order parameter shows more 
inhomogeneous behavior with larger kpS^Q. 



fa}: 



W'Fj 5 




FIG. 1. (Color online) Contour plots of the order parameter 
A(a;,y) with the corresponding diagonal profiles in the ab- 
sence of applied magnetic field. Panels (a)-(c) are for sample 
I with sizes W/^o ~ 5, 8, 10, respectively. Panels (d) and (e) 
are for sample II with sizes 5^o and lO^o, respectively. Panel 
(f) corresponds to sample III with size 5^o- 




FIG. 2. (Color online) Size dependence of the order parameter 
in the center of the sample |Ac| and the spatially averaged 
I A ml for sample I. 
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FIG. 3. (Color online) Free energy as a function of the magnetic flux through sample I with size lO^o x lO^o. Different colors 
indicate the different winding numbers L and the shaded area indicates the flux range over which the vortex state with winding 
number L is the ground state. The vortex configurations of the ground state marked by solid squares are shown separately in 
Fig. |4] and the corresponding free energy curves are marked by dots. The inset zooms in the region where L = h and 6 are the 
ground states 



IV. VORTEX CONFIGURATIONS AT ZERO 
TEMPERATURE 

In this section, we consider the zero temperature case 
for which the system is always in the quantum limit 
(T < l/fci?^o)- First, we study the sample I with size 
10^0 X 10^0 and show in Fig.[3]the free energy of the sta- 
ble vortex states for flux $/<i>o G [0, 20]. Different curves 
(colors) represent states with different winding number 
L and the states among them which reached the ground 
state are marked by dots. Vertical lines and shadows 
show the flux range for each L state as the ground state. 
The top dashed line stands for the free energy of the sys- 
tem in the normal state when the coupling constant g is 
set to zero. When compared with the GL theory^^, one 
of the differences is that the free energy of the normal 
state depends on the magnetic field while it is a constant 
in GL theory. The reason is that the energy levels of 
the confined electrons are different for different magnetic 
fields. In our case, the change of the energy is relatively 
small when compared with the energy gap (energy dif- 
ference between the normal state and the superconduct- 
ing state) especially in weak fields. Although the shown 
energy curves look conventional, there are significant dif- 
ferences with the GL theoryiii 

By sweeping the magnetic field up and down, we can 
get the full energy spectrum and the corresponding vor- 
tex states. For a certain magnetic field, it is common to 
have more than one converged solution. The lowest en- 
ergy state is the ground state while the states with higher 
energy are referred to as metastable states. In Fig. HI we 



show the contour plot of the absolute value of the order 
parameter of the corresponding ground states for various 
winding numbers. 

As can be seen from Fig. [31 the system favors states 
with winding numbers L = 1,4,8 and 9 because they 
have relatively large ground state flux range (excluding 
the Meissner state). From Fig. |4l we observe that these 
states have fourfold symmetry which is compatible with 
the sample geometry. One interesting feature of this sys- 
tem is the richness of metastable states. These states ap- 
pear for all winding numbers L except L = and 1 . The 
number of metastable states reaches a peak for L = 4 
and equals 11. From the free energy curves, we notice 
that the energy difference between the ground state and 
the metastable states can sometimes be very small. This 
makes the ground state difficult to find in simulations 
unless we sweep the field up and down many times. We 
also note that most metastable states are focused at lower 
magnetic fields from the corresponding ground state flux 
range. The reason is that vortices get easily stuck at the 
boundary due to the pronounced oscillations of the order 
parameter. For the same reason, their stability range is 
narrower due to asymmetry. 

The number of metastable states decreases for higher 
L. In this case, the shorter distances between vortices 
cause strong interaction between them. This limits the 
choice of stable positions for vortices and therefore the 
number of metastable states are lower. Due to this rea- 
son, metastable states are less favorable for smaller sam- 
ples because of easy saturation with vortices. For exam- 
ple, no metastable states were found in sample I with size 
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FIG. 4. (Color online) Contour plot of the amplitude of the 
order parameter for the ground states marked by solid squares 
shown in Fig. [3l L = 10(q:) and are for flux $/$o = 18.2 
and 18.7, respectively. Darkest areas indicate the positions of 
vortices. Note that at the center an anti-vortex (for L = 7) or 
a giant vortex [for L — 10(/3)] can spontaneously form. Note 
also that for some vorticities two different ground states are 
found (at different magnetic fields). 
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Ground states 



Next we discuss the ground states configurations for 
different vorticities in more detail. For the L = 1 state, 
shown in Fig. 21 the vortex sits at the center of the sample 
which is compatible with the conventional picture. Al- 
though such result is observed for different parameters of 
the sample, the state with diagonal location of the vortex 
can also be found in some cases.^"^ It is clear that the or- 
der parameter around the vortex core is suppressed and 
the profile shows the competition between C4 symmetry 
and Coo ■ It means that the vortex has long range (longer 
than ^0) interaction with other objects such as the other 
vortices and/or boundaries. 
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FIG. 5. (Color online) Phase diagram for samples I-III. Dif- 
ferent shadings of the background indicate different types of 
vortex states. Black(blue) dots indicate vortices, and dia- 
monds indicate antivortices in the schematic diagrams of vor- 
tex configurations (bottom figures). The symbols are larger 
when vortices contain multiple flux quanta. 
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By increasing the flux to $/$o = 5.24, the ground 
state shifts from L = 1 to L = 2 and two vortices sit along 
the diagonal. This again coincides with the result from 
GL theory, and results from the competition between the 
confinement imposed on vortices by the Meissner currents 
and the vortex- vortex repulsion. In this case, these effects 
can be clearly seen from the profile of the order parame- 
ter, especially from the suppressed area around the vor- 
tices. The vortex- vortex interaction suppresses the order 
parameter mostly in the area between them. This can 
not be found in GL theory because the order parameter 
is always smooth and changes slowly in space. The state 
with L = 3 shown in Fig. |4] becomes the ground state 
when the fieldapplied flux is between $/$o = 6.9 and 
8.12. It resembles the multi- vortex state obtained within 
GL theory where the three vortices are at the apices of a 
equilateral triangle. However, the perpendicular bisector 
of the triangle always coincides with one of the diagonals 
of the square sample in GL theory while it is parallel to 
one of the edges in our case. This is because the grid- 
like pattern in the inhomogeneity of the order parameter 
imposes preferential positions for the vortices inside the 
square. The state L = 4 has a similar feature but the 
configuration is compatible with the GL result. 

From these states, we conclude that when the GL vor- 
tex configuration, which minimizes the vortex- vortex and 
vortex-boundary interaction, matches the oscillation pat- 
tern due to quantum confinement, then the state has a 
wider flux stability range. 

Two ground states, L = h{a) and L — 5(/3), are found 
for L = 5 in the flux ranges $/$o = 10.48 - 10.78 and 
10.78 — 11.62, respectively. Both of them have a pentag- 
onal vortex configuration. This is because the particular 
shape resonance at the considered field causes the order 
parameter to be peaked at the center. Therefore, it costs 
energy for vortices to sit in the center of the sample. 
For the same reason, the ground states L = 6(a) and 
L = 6(/3) do not have vortices in the center. Moreover, 
when L ^ 6, vortices start to be compressed in the sam- 
ple. If they do not form giant vortices, they will be very 
close to each other and form string-like structures [see 
L = 6(a) state]. 

States with L = 7, L = 8, L = 9 and L = 10(;3) shown 
in the panels of Fig. |4] have a common feature, as all 
of them keep the fourfold symmetry. L — 7 contains a 
antivortex at the center while L = 9 and L — 10(/3) have 
a single vortex and a giant vortex with 2<l'o a-t the center, 
respectively. The state with L = 7 is the only ground 
state which contains an antivortex. The antivortex is 
closely surrounded by four vortices and forms the core 
structure for L — 7. The outer shell is formed by the 
remaining four individual vortices sitting at four corners. 
States with L — 8 and L = 9 contain vortex dimers, 
i.e. two vortices close to each other at each corner. The 
fourfold symmetry makes both former states have a larger 
ground state flux range. L — 10(/3) also keeps the C4 
symmetry but the energy is sometimes even higher than 
the state L = 10(a), which has only C2 symmetry. The 



reason is that the giant vortex costs extra energy. 

In order to visualize the changes in the ground states 
when key parameters change, we plot the phase diagram 
for samples I-III for W/£,o = 5 - 10 and $/$o = - 10 
in Fig. O Different shadings of blocks in Fig. [S] indicate 
different vortex types. The plain white background rep- 
resents conventional multi-vortex states as found within 
the GL theory, while the blue background with square 
grid represents giant vortices, also compatible with the 
result obtained from GL theory. Asymmetric vortex 
states attained only by BdG theory are represented by 
yellow background with horizontal grid pattern. States 
containing parallel vortex chains, represented by orange 
background with vertical grid pattern, and part of the 
vortex-antivortex molecules represented by pink(grey) 
background, are new compared to GL theory. 

We can conclude from Fig. [5] that the quantum size ef- 
fect is important not only for the transition field and the 
stability range of the different vortex states, but also for 
the vortex configurations. The reason is that the oscilla- 
tion patterns of the order parameter are very sensitive to 
kp^o and may cause totally different behaviors even for 
two samples of identical size. 

As can be seen from Fig. [H^a) and (b) , the vortex phase 
transition fields vary greatly with sample size W except 
for the phase boundary between L = and L = 1 state. 
All the phase boundaries oscillate with W. We find that 
some samples favor vortex states with even winding num- 
ber L while other disfavor them. For example, sample I 
with W/S^o — 5 favors L — 2 state whereas the one with 
W^/ £,0 ~ T disfavors L — 4 to the point of non-existance. 
When W becomes large, the transition fields start to con- 
verge and the flux stability range of each L state will be 
roughly one flux quantum. 

When compared to samples I and II, sample III (with 
large kpS^o) shows a more conventional picture. More- 
over, the phase transition held increases only slightly 
with increasing W. It means that for large kp^o, the 
quantum size effect on the transition field, at least when 
winding number L is small, can be neglected. Never- 
theless, a plethora of different vortex configurations is 
found. For example, an asymmetric L ~ 1 state is found 
in all three samples. Moreover, sample I always shows 
asymmetric states when W/£,o < 6. One other interest- 
ing phenomenon is that, for L — 2, sample II can host 
all five types of vortex states with W increasing. 

From the phase diagram in Fig. [5l we notice that 
nanoscale superconductors favor anti-vortices and disfa- 
vor giant vortices. For example, the giant vortex state 
appears for L = 2 in sample II with W/£,o = 7. Based 
on GL theory, only smaller samples will exhibit a giant 
vortex configuration. However, when W/£q < 7, the two 
individual vortices form an asymmetric vortex state. On 
the other hand, we find that the probability of forming 
anti-vortices is much higher than in GL theory. In GL 
case, antivortex states usually appear for L = 3 when 
four vortices are at the four corners and surround a cen- 
tered antivortex. Usually the distance between vortices 
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FIG. 6. (Color online) Free energy curves of states with 
L = 2 in sample I with size W/£,o ~ 10 (top panel) and 
the corresponding contour plots of the order parameter (bot- 
tom panel). Note that the contour plot of the ground state 
L = 2(e) is shown in Fig.[4l The inset shows the energy differ- 
ence between the ground state L — 2(e) and the metastable 
states L = 2(f) and L — 2(g) at higher field. 



and the antivortex is small (less than ^). In the BdG 
calculation, at least two more antivortex states can be 
found. One is still the L = 3 state and the antivortex 
is still at the center but four vortices are at the edges 
instead of the corners. This configuration was first pre- 
sented by us in Ref. js^ where we showed that the size 
of the vortex/anti-vortex(V-aV) molecule is larger than 
the one obtained with GL theory. The other antivortex 
state appears for L = 2. In this case, the four vortices 
still sit at four corners but the centered antivortex car- 
ries two flux quanta, e.g. it is a giant antivortex! Due to 
the strong vortex-antivortex interaction, the size of such 
V-aV molecule is small. 



B. Metastable states 

Next we will briefly discuss the metastable states of 
this system. To do this, we start from sample I with 
W/£,o = 10. Metastable states are important in the 
BdG formalism because the energy difference between 
the ground and the metastable states can be very small. 
This suggests that these states could be easily found in 
experiments. Alternatively, some metastable states can 
become ground states as the parameters are changed. 

All six found metastable states for L = 2 and their free 
energy curves are shown in Fig. [Bl The state L — 2(/) 
is similar to the ground state L = 2(e), but rotated over 
45°, hence their free energies are very close to each other. 
Actually, the difference in the orientation of the vortex 
pattern always results in a small difference in energy. 
State (f) is not obtained within the GL theory. In our 
case, due to the shape-resonant inhomogeneity of the or- 
der parameter, the rotation of the vortex pattern to the 
ground-state configuration is prevented by the spatial os- 
cillations of the order parameter. 

The metastable state L = 2{g) is only stable at higher 
field and its free energy is very close to the ground state 
L = 2(e). Therefore we zoomed on the energy difference 
in the inset of Fig. |6l From the figure, one sees that the 
energy of the L — 2{g) state is lower than the ground 
state, L — 2(e), when the applied fiux is larger than 
^/^Q = 7.3. In fact, the state can exist even up to 
^/^o — 10. From the vortex configuration shown in 
Fig. [HI we find that this state is a giant vortex. Such a 
state has been predicted by the GL theory because the 
magnetic field pushes the two vortices towards each other 
and makes them merge into a giant vortex. Usually, the 
phase transition between the multi vortex state and the 
giant vortex state is continuous (second-order). However, 
the barrier induced by the inhomogeneity of the order 
parameter leads to a first order phase transition in our 
case. One more difference between the BdG giant vortex 
and the one in GL theory is its core structure. Due to 
the shape resonances, the contour plot of the core shows 
a diagonal cross shape while the giant vortex core in the 
GL case is always circular. Furthermore, the giant vortex 
state in our results has two allotropes: see state L = 
2(c), compared to the state L = 2{g). The L = 2(g) 
state exists up to higher field while L = 2(c) only exists 
in lower field. Hence, the size of the giant vortex seen 
in L = 2(c) is larger than the one seen in L = 2{g). 
Another difference between them is the orientation of the 
core. L = 2{g) has diagonal cross shape while the state 
L = 2(c) has edge cross shape. 

Other three metastable states L = 2(a), L = 2(6) and 
L — 2(d) are observed only in lower field. They have in 
common the fact that at least one vortex is stuck at the 
boundary since the Meissner current pushes the vortex 
outward at low fields. It is obvious that they have lower 
energy when the applied field is lower. The energy of 
state L = 2(a) is always lower than the one of state L = 
2(6) because of the longer distance between the vortices. 
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FIG. 7. (Color online) Free energy curves of states with L — 
3 in sample I with size W/^o = 10 and the corresponding 
contour plots of the order parameter. The contour plot of the 
ground state L = 3{g) is shown in Fig.H 



As is usual in mesoscopic superconductors, vortices in 
these states avoid to be located at the very corners of 
the sample, due to strong local superconductivity there. 

For the L — 3 metastable states, the results are sum- 
marized in Fig. [71 The state L = 3{h) is a V-aV state and 
exists only in higher field while the giant vortex L = 3(c) 
only exists in lower field. Note again that the energy 
of the giant vortex states is much higher than the other 
metastable states with three single vortices. L = 3(/) has 
lowest energy for L = 3 around $/$o = 6, when there 
is one vortex located at the boundary. States L = 3(6) 
and L = 3(a) follow when the field decreases and there 
are two and three vortices stuck at the boundaries, re- 
spectively. States L — 3(d) and L = 3(e) are disfavored 
and have higher energy due to the close distance between 
vortices. Note again that no vortex sits at the corners in 
this states. 

States with L — 4 show a wide ground state flux range 
and 11 different metastable states, which is the largest 
variety of all L vortex states. From Fig. [SI we find that 
the metastable states concentrate around the applied fiux 
$/$o = 6.5. 
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FIG. 8. (Color online) Free energy curves of vortex states 
with L — 4 in sample I with size W/S,o — 10 and the contour 
plots of the order parameter of selected vortex states. Inset 
shows details of the free energy curve around $/$o = 6.5. 
The contour plot of ground state L — 4{k) is shown in Fig.H 



For low fields, we conclude again that vortices are close 
to the surface and these states are always the lowest 
energy state for a given L state at low fields. States 
L = 4(5) and L = 4(c) have C4 symmetry and all of the 
four vortices are trapped close to the boundary. Please 
note that in the state L = 4(a) four vortices sit at the 
corners. This kind of state is rare in the BdG results 
because corners give the highest potential energy contri- 
bution for vortices. From the free energy curve of the 
state L = 4(a), we can find that the slope of the en- 
ergy curve is opposite to the other L — 4 states [such as 
L = 4(5)] in this field range. This indicates that vortices 
are repelled by the Meissner current in order to balance 
the inward force the vortices experience from the corner. 
When the field is too low, a vortex is expelled from the 
sample and the state jumps to a L = 3 state. The vortex 
configuration L — 4{f) has been found experimentally 
in conventional mesoscopic superconductors* but was a 
result of the presence of pinning sites. This state can not 
be obtained in plane squares within the GL theory . An- 
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other state to notice is L = 4(z) whose energy curve does 
not cross any other L — 4: curve. When the supercon- 
ductor is in this state and the field is swept down, this 
state will be the first to jump to the L — 3 state. This is 
understandable from the vortex configuration of L = 4(i) 
because the vortex at the corner is easily expelled when 
field is lowered. 

At high fields only one metastable state exists L = 
4(1). It can be seen as the state obtained after a 90 
degrees rotation of the ground state-L = 4(A:). This is a 
consequence of the fact that the inhomogeneous pattern 
of the order parameter changes with field. At such a 
high field, the corner vortex position in state L = A{k) 
becomes unstable, which forces the vortices to sit at the 
edges, similar to the L = 4{l) case. At the same time, 
the strong field pushes vortices closer together so that 
the distance between vortices in L = 4(Z) is shorter than 
the one in the ground state L — 4(fc). 



V. VORTEX STATES AT FINITE 
TEMPERATURE 

So far, all our calculations were done at zero tempera- 
ture, T/Tc = 0, where Tc is the bulk critical temperature 
at zero fiux $/$o = 0. In what follows, we investigate 
the effect of temperature on the vortex configuration. 
First, we show all the vortex states for the fiux range 
<f>/<f>o e [0,10], for sample I with size W/S,o = 10 at 
T/Tc = 0.6 where the system is NOT in the quantum 
limit since T/Tc > l/^i="Co- The corresponding free en- 
ergy curve as a function of flux is presented in Fig. 1^1 
Contrary to the results for T/Tc — 0, which were shown 
in Fig. [3l the figure looks more conventional (similar to 
the results obtained by GL theory in Ref . TT) and there is 
only one stable state for each winding number L. More- 
over, only giant vortex states are found in this case for 
L ^ 2. For the size of the square sample considered 
here, W « 10^, GL theorjtii predicts that multi-vortex 
states should exist. Here we find instead that multi- 
vortex states are absent since ^ increases as temperature 
increases. 

In order to see how temperature affects the coherence 
and the profile of the superconducting order parameter, 
we show the order parameter for <I>/<I>o = 4 and L = 1 in 
Fig. [ini The diagonal profile of |A| at T/Tc = shows 
the strongest Friedel-like oscillations. As temperature 
increases to T/Tc ~ 0.2, the profile is similar to the one 
obtained at zero temperature, but with less oscillations at 
the vortex core. Both cases are in the quantum limit and 
the order parameter shows rapid variation in the core. 
When temperature reaches T/Tc = 0.6, we find that both 
the average and the oscillations of the absolute value of 
the order parameter are suppressed which indicates that 
the vortex states become more conventional. Finally, the 
order parameter is smooth at T/Tc = 0.8 and the GL 
results are approached. 

As seen from the Fig. [TOl the coherence length, which 




FIG. 9. (Color online) Free energy as a function of the mag- 
netic flux through the square sample I for size W/^o ~ 10 and 
T/Tc = 0.6. 

represents the vortex core radius, increases with increas- 
ing temperature. As defined by Kramer and Pesch^, we 
calculate the coherence length as 



where r is the distance to the vortex core. We plot in 
Fig. [TT] ($i/$o)~^ as a function of temperature, T/Tc ■ 
As discussed in Ref. [s^, can be described by ^(T) oc 
{Tc - T)-i/2 when T is close to Tc {T/Tc > 0.5 in our 
case) . In the intermediate temperature regime, there is a 
substantial suppression of the coherence length because 
of the bound states. At low temperature, the shrinkage 
of the coherence length stops and saturates when the 
system is in the quantum limit. Note that at T/Tc = 
0.6 is around three times larger than the one at zero 
temperature. This explains why only giant vortex states 
can be found at such temperatures. 

Fig. [T^] shows the order parameter for sample I at 
T/Tc = 0.6 for L = 2 in panel (a) and L = 3 in panel 
(b), respectively. Both are giant vortex states and the C4 
symmetry grid pattern is strongly suppressed. As can be 
seen from the figure, the vortex cores show perfect cir- 
cular symmetry, which is in agreement with the results 
from GL theory. Of course, the size of the vortex core 
shown in panel (b) is larger than the one shown in panel 
(a) because its vorticity is larger. 

Finally, we end this section with the (T-$) phase dia- 
gram for lower fields for sample I with W/^q — 5. This 
is shown in Fig. [T31 The thick black curve indicates 
the phase boundary between the superconducting and 
the normal state. When the system is in the quantum 
limit, for these parameters, only unconventional vortex 
states, such as asymmetric L = 1 and L = 3 states and 
edge-parallel L = 2 states, are found as ground states. 
When temperature increases, the vortex states become 
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FIG. 10. (Color online) The diagonal profile of the order 
parameter for sample I with size W/£,o = 10 for $/<i>o ~ 4 and 
L = 1 at different temperatures. The corresponding contour 
plots of the order parameter at T/Tc = and T/Tc = 0.6 are 
also shown. 
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FIG. 12. (Color online) The contour plots of the order pa- 
rameter for sample I with size W/^o = 10 at T/Tc = 0.6 for 
(a) L = 2 and (b) L = 3. 



Concluding this section, higher temperature: 1) makes 
vortex states look more conventional (closer to the GL 
results); 2) smoothens the order parameter; 3) suppresses 
the influence of the oscillation of the order parameter and 
4) increases the superconducting coherence length ^. As 
a consequence, the number of metastable states is also 
lowered. The effect of temperature is very different (more 
complex) from the effect obtained by simply changing the 
effective size of the sample as is usually done within the 
GL theory. 



VI. GIANT ANTI- VORTEX AND THE 
STRUCTURE OF THE VORTEX CORE 

In this section, we discuss the appearance and stability 
of anti- vortex states in the BdG theory in order to explain 
the existence of the giant antivortex. Actually, such a 



FIG. 11. (Color online) Temperature dependence of 
l/(Ci/Co)^ for sample I with size VK/^o = 10 and $/<l>o = 4. 

conventional and the C4 symmetry of the states is always 
preserved. Note that the asymmetric L = 1 state goes 
through a continuous phase transition to the symmetric 
L = 1 state, which means the vortex moves gradually as 
the temperature changes. However, for higher winding 
numbers, the system usually goes through a first order 
phase transition. For example, the phase transition be- 
tween the parallel vortex state and the giant vortex state 
of L = 2 is of first order. This is different from the GL 
result, where vortices merge into a giant vortex through 
a continuous phase transition^. For the L = 3 state, 
we note that the ground state flux range for the four- 
fold symmetric V-AV state is larger than the asymmetric 
one due to the compatibility of its symmetry with the 
geometry of the sample. 
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FIG. 13. (Color online) Temperature-Flux phase diagram for 
sample I with size W/^o = 5. The vortex configurations of 
areas (i-iii) are shown as insets in the upper right corner. 
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FIG. 14. (Color online) Diagonal profile of the order param- 
eter for different vortex states in the GL theory. Blue thick 
curves represent the order parameter. All the phases have 
been adjusted such that on the diagonal the order parameter 
is real. The horizontal line indicates the zero of the order 
parameter and a vortex appears when it intersects the order 
parameter. Panels (a) and (c) represent multi-vortex states. 
Panels (b) and (d) are giant vortex states and (e) is a vortex- 
antivortex configuration. 



state was already found in Ref.^^ through the hnear GL 
method by introducing artificial pinning. 

From the phase diagram shown in Fig.[5l we found that 
anti- vortex states are surprisingly stable within BdG the- 
ory. This is due to the fact that the grid pattern oscil- 
lation of the order parameter gives an additional contri- 
bution to the symmetry of the vortex states and there- 
fore, in a square sample, the C4 symmetry is enhanced. 
The other reason to form an anti-vortex is that the os- 
cillations induced by the order parameter are seen in the 
vortex core where the order parameter is already sup- 
pressed. These oscillations can easily lead to a shift in 
the phase of the order parameter by tt and, thus, result 
in the formation of vortex-antivortex molecules. 

In order to explain this, we first discuss briefly the vor- 
tex profile in the GL theory. Fig. [T3] shows a schematic 
diagram of the vortex states for different winding number 
L in GL theory. The diagonal profiles of the order pa- 
rameters vary smoothly in space and the vortex emerges 
where the order parameter vanishes. Note that the phase 
of the order parameter is adjusted such that along the di- 
agonal the order parameter is real. Panel (a) from Fig. [HI 
shows the simplest case when only one vortex sits at the 
center. As can been seen from the figure, the order pa- 
rameter changes sign, which indicates the tt phase shift 
of the order parameter. The profile is an odd function 
and A(r) ^ r near the vortex core. Panels (b) and (c) 
from Fig. [HI show the diagonal profiles for L = 2. Both 
profiles are even functions due to the 27r phase shift be- 
tween the opposite corners. The order parameter exhibits 



FIG. 15. (Color online) Similar as Fig. [14] but now for the 
BdG theory. Panels (a) and (b) show the symmetric and 
asymmetric L = 1 vortex states, respectively. Panels (c) and 
(d) are the giant vortex and vortex-giant antivortex L = 2 
states, respectively. Panels (e) and (f) are the giant vortex 
and vortex-antivortex L — 3 configurations, respectively. 




A(r) 



property. When there is only one root, as can 



be seen from panel (b), the vortex is a giant one. When 



FIG. 16. (Color online) Vortex state for sample III with 
W/(,o = 7 at ^/^o = 20 with winding number L = 12. Panels 
(a) and (b) show contour plots of the order parameter and its 
phase, respectively. Panel (c) shows schematically the vortex 
configuration. Dark(blue) dots and open diamonds indicate 
vortices and antivortices, respectively. 



there are two roots, as shown in panel (c), the configura- 
tions are multi- vortex states. Similarly, the profiles of the 
order parameter shown in panels (d) and (e) from Fig. [14] 
for L = 3 show a A(r) ^ spatial dependence. One root 
means that we have a giant vortex state whereas three 
roots represent a vortex-antivortex configuration. Note 
that, in order to generate the central anti- vortex, the or- 
der parameter has to oscillate around the center of the 
square. 

Now let us move to nano-size superconductors where 
the BdG theory has to be used and the spatial oscillation 
of the order parameter cannot be neglected. As can be 
seen from Fig. [15] the oscillation plays an important role 
in generating vortices, especially when the value of the 
order parameter is comparable to the amplitude of the 
oscillation. For instance, panels (a) and (b) from Fig. [15] 
show the symmetric and the asymmetric L — 1 vortex 
states. The reason for the appearance of the asymmetric 
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vortex state is the fact that the order parameter has an 
odd number of osciUations across the diagonal Thus, the 
vortex can not sit at the center. For the L = 2 states, 
panel (c) shows a giant vortex state where the sign of the 
profile of the order parameter is always positive across 
the diagonal. However, due to the oscillations, the case 
shown in panel (d) of Fig. [15] can easily exist and shows 
a giant anti-vortex (L = —2) at the center. Further, the 
configurations show a large diversity for a fixed winding 
number L. When L = 2, the configuration can be 1 + 1, 
2 + 0, 3 — 1, 4 — 2 and so on. Panels (e) and (f) from 
Fig.[T5]are for L ~ 3 states. Apparently, they are similar 
to the GL case shown in Fig. [141 but the probability of the 
occurrence of the V-aV state is much larger than in GL 
case. The reason is that the result with one root is just 
a special case while the general case shows oscillations at 
and around the vortex core. 

The V-aV molecules do not only exist for smaller wind- 
ing number L, but they can also appear for large L in 
the BdG results. Fig. \W\ shows an example for sam- 
ple III with W^/^o=7 at <i>/$o = 20 and with a winding 
number L = 12. We find that vortices concentrate in 
the central dark (blue) area where the order parameter 
is strongly suppressed. From the phase of the order pa- 
rameter, which is shown in Fig. I16f b). the total winding 
number L = 12 is found but it is difficult to distinguish 
each vortex. After a careful analysis, we plot schematic 
diagram of the vortex configuration in Fig. I16f c). The 
dark(blue) dots and open diamonds indicate vortices and 
anti- vortices, respectively. As can be seen, the lattices of 
vortices and antivortices are nested within each other. 
Since anti-vortices attract vortices, all the vortices (24 
vortices and 12 antivortices) can be condensed in the cen- 
tal area of the sample. This picture becomes more accu- 
rate when kp^Q is large. The stronger the oscillations of 
the order parameter the more V-aV pairs are generated. 
However, the size of the V-aV pair can only be of the or- 
der of the Fermi wavelength. Thus, it will be very hard 
to detect them in experiments. This is why these states 
are mostly treated as a giant vortex in conventional su- 
perconductors. In other words, the suppressed central 
area of the order parameter, after coarse graining, will 
look like a giant vortex with L = 12. 



VII. CONCLUSION 

To summarize, we investigated the vortex states in a 
nanoscale superconducting square for difi'erent sizes W, 
parameters fcF^o, and temperatures T. First, we found 
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that the inhomogeneous pattern of the order parameter 
in the absence of magnetic field strongly depends on kp£^Q 
and the size W. This oscillation pattern will give an ad- 
ditional contribution to competing effects that determine 
the vortex configurations when the field is applied. Due 
to the inhomogeneous order parameter induced by the 
quantum topological confinement, samples with different 
kpC^o and W will favor different winding numbers L. 

We find unconventional vortex states such as asym- 
metric, edge-parallel and vortex-antivortex states as the 
ground state of our nanoscale system. These were never 
seen in the Ginzburg-Landau approach. The inhomo- 
geneous pattern of the order parameter, especially the 
strong oscillation at the boundaries causes additional po- 
tential wells for vortices which in turn generates a lot of 
metastable vortex states. Furthermore, in the quantum 
limit, nano-size superconductors favor vortex-antivortex 
molecules while disfavoring giant vortex states. 

We observe that vortex ground states and the phase 
transition fields are very sensitive to changes in the pa- 
rameter fc_F5oi size W and temperature T. This is a direct 
consequence of the quantum size effect. However, this 
effect is suppressed when the size W is large or when 
temperature is high. In this case most metastable states 
become unstable and the ground states become compat- 
ible with GL theory. 

For high magnetic fields, vortex-antivortex pairs can 
be easily found when kp^o is large because the absolute 
value of the order parameter becomes smaller than the 
amplitude of its oscillations. However, detection of such 
states is beyond the current experimental abilities. 

The peculiar vortex states uncovered in the present 
work should be observable in superconducting systems 
where kp^^Q is small. Such systems could be high-Tc 
superconducting nano-grains for which the coherence 
length is small or cold-atom condensates with small kp, 
i.e. large Fermi wavelength. Of special interest could 
be hybrid systems made of superconducting substrates 
and graphene sheets for which the Fermi wavelength is 
highly tunable near the Dirac point. Future work could 
also address the fundamental vortex-vortex and vortex- 
antivortex interactions for systems with a small kp^o, for 
which the oscillations of the order parameter on the order 
of Aj? become important. 
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